Scedastic Surrogate Swarm R Package

rm(list=ls())
library(ScedasticSurrogateSwarmR)
library(keras)
library(dplyr)
set.seed(128)

Problem Setup

In the previous example, SSSoR was used to optimize the weights of a neural network. In that example the simple neural network was built from scratch with matrix operations, but it is important to note that SSSoR can also optimize models built in Keras. This also gives us the opportunity to compare performance to a gradient-based optimizer like Adam.

ndat = 100
xdat=sort(runif(ndat,-10,10))
ydat = c(sin(0.5*xdat)/2)+0.6

Optimization of Keras Model with Adam

We will use keras to build an identical neural network to our previous example.

# Define the model
model <- keras_model_sequential()

# Add a single hidden layer with 64 neurons and ReLU activation
model %>%
  layer_dense(units = 20, input_shape = c(1)) %>%
#  layer_dense(units = 20, input_shape = c(1)) %>%
  layer_activation_leaky_relu(alpha = 0.03)

# Add the output layer with a single neuron (regression problem)
model %>%
  layer_dense(units = 1,use_bias = F)

# Compile the model
model %>% compile(
  loss = 'mean_squared_error', # You can choose a different loss function if needed
  optimizer = optimizer_adam() # You can choose a different optimizer if needed
)

We will use Adam to train this neural network.

# Fit the model to your data
history <- model %>%
  fit(
    x = xdat,
    y = ydat,
    epochs = 100, # You can adjust the number of epochs as needed
    batch_size = 1,
    verbose = 2
  )
#> Epoch 1/100
#> 100/100 - 1s - loss: 6.3931 - 1s/epoch - 12ms/step
#> Epoch 2/100
#> 100/100 - 0s - loss: 1.8310 - 318ms/epoch - 3ms/step
#> Epoch 3/100
#> 100/100 - 0s - loss: 0.4983 - 278ms/epoch - 3ms/step
#> Epoch 4/100
#> 100/100 - 0s - loss: 0.2165 - 315ms/epoch - 3ms/step
#> Epoch 5/100
#> 100/100 - 0s - loss: 0.1762 - 262ms/epoch - 3ms/step
#> Epoch 6/100
#> 100/100 - 0s - loss: 0.1596 - 267ms/epoch - 3ms/step
#> Epoch 7/100
#> 100/100 - 0s - loss: 0.1507 - 237ms/epoch - 2ms/step
#> Epoch 8/100
#> 100/100 - 0s - loss: 0.1383 - 326ms/epoch - 3ms/step
#> Epoch 9/100
#> 100/100 - 0s - loss: 0.1285 - 218ms/epoch - 2ms/step
#> Epoch 10/100
#> 100/100 - 0s - loss: 0.1187 - 234ms/epoch - 2ms/step
#> Epoch 11/100
#> 100/100 - 0s - loss: 0.1062 - 224ms/epoch - 2ms/step
#> Epoch 12/100
#> 100/100 - 0s - loss: 0.0949 - 218ms/epoch - 2ms/step
#> Epoch 13/100
#> 100/100 - 0s - loss: 0.0875 - 229ms/epoch - 2ms/step
#> Epoch 14/100
#> 100/100 - 0s - loss: 0.0774 - 219ms/epoch - 2ms/step
#> Epoch 15/100
#> 100/100 - 0s - loss: 0.0722 - 218ms/epoch - 2ms/step
#> Epoch 16/100
#> 100/100 - 0s - loss: 0.0660 - 229ms/epoch - 2ms/step
#> Epoch 17/100
#> 100/100 - 0s - loss: 0.0594 - 232ms/epoch - 2ms/step
#> Epoch 18/100
#> 100/100 - 0s - loss: 0.0563 - 217ms/epoch - 2ms/step
#> Epoch 19/100
#> 100/100 - 0s - loss: 0.0473 - 217ms/epoch - 2ms/step
#> Epoch 20/100
#> 100/100 - 0s - loss: 0.0411 - 218ms/epoch - 2ms/step
#> Epoch 21/100
#> 100/100 - 0s - loss: 0.0368 - 232ms/epoch - 2ms/step
#> Epoch 22/100
#> 100/100 - 0s - loss: 0.0368 - 232ms/epoch - 2ms/step
#> Epoch 23/100
#> 100/100 - 0s - loss: 0.0323 - 218ms/epoch - 2ms/step
#> Epoch 24/100
#> 100/100 - 0s - loss: 0.0266 - 233ms/epoch - 2ms/step
#> Epoch 25/100
#> 100/100 - 0s - loss: 0.0253 - 231ms/epoch - 2ms/step
#> Epoch 26/100
#> 100/100 - 0s - loss: 0.0220 - 232ms/epoch - 2ms/step
#> Epoch 27/100
#> 100/100 - 0s - loss: 0.0206 - 216ms/epoch - 2ms/step
#> Epoch 28/100
#> 100/100 - 0s - loss: 0.0177 - 218ms/epoch - 2ms/step
#> Epoch 29/100
#> 100/100 - 0s - loss: 0.0174 - 219ms/epoch - 2ms/step
#> Epoch 30/100
#> 100/100 - 0s - loss: 0.0158 - 235ms/epoch - 2ms/step
#> Epoch 31/100
#> 100/100 - 0s - loss: 0.0144 - 232ms/epoch - 2ms/step
#> Epoch 32/100
#> 100/100 - 0s - loss: 0.0132 - 234ms/epoch - 2ms/step
#> Epoch 33/100
#> 100/100 - 0s - loss: 0.0116 - 218ms/epoch - 2ms/step
#> Epoch 34/100
#> 100/100 - 0s - loss: 0.0114 - 232ms/epoch - 2ms/step
#> Epoch 35/100
#> 100/100 - 0s - loss: 0.0118 - 218ms/epoch - 2ms/step
#> Epoch 36/100
#> 100/100 - 0s - loss: 0.0100 - 247ms/epoch - 2ms/step
#> Epoch 37/100
#> 100/100 - 0s - loss: 0.0099 - 218ms/epoch - 2ms/step
#> Epoch 38/100
#> 100/100 - 0s - loss: 0.0098 - 216ms/epoch - 2ms/step
#> Epoch 39/100
#> 100/100 - 0s - loss: 0.0097 - 230ms/epoch - 2ms/step
#> Epoch 40/100
#> 100/100 - 0s - loss: 0.0082 - 235ms/epoch - 2ms/step
#> Epoch 41/100
#> 100/100 - 0s - loss: 0.0085 - 219ms/epoch - 2ms/step
#> Epoch 42/100
#> 100/100 - 0s - loss: 0.0080 - 219ms/epoch - 2ms/step
#> Epoch 43/100
#> 100/100 - 0s - loss: 0.0075 - 231ms/epoch - 2ms/step
#> Epoch 44/100
#> 100/100 - 0s - loss: 0.0076 - 220ms/epoch - 2ms/step
#> Epoch 45/100
#> 100/100 - 0s - loss: 0.0068 - 229ms/epoch - 2ms/step
#> Epoch 46/100
#> 100/100 - 0s - loss: 0.0078 - 217ms/epoch - 2ms/step
#> Epoch 47/100
#> 100/100 - 0s - loss: 0.0067 - 237ms/epoch - 2ms/step
#> Epoch 48/100
#> 100/100 - 0s - loss: 0.0067 - 251ms/epoch - 3ms/step
#> Epoch 49/100
#> 100/100 - 0s - loss: 0.0072 - 223ms/epoch - 2ms/step
#> Epoch 50/100
#> 100/100 - 0s - loss: 0.0065 - 231ms/epoch - 2ms/step
#> Epoch 51/100
#> 100/100 - 0s - loss: 0.0058 - 218ms/epoch - 2ms/step
#> Epoch 52/100
#> 100/100 - 0s - loss: 0.0062 - 234ms/epoch - 2ms/step
#> Epoch 53/100
#> 100/100 - 0s - loss: 0.0061 - 235ms/epoch - 2ms/step
#> Epoch 54/100
#> 100/100 - 0s - loss: 0.0055 - 217ms/epoch - 2ms/step
#> Epoch 55/100
#> 100/100 - 0s - loss: 0.0056 - 218ms/epoch - 2ms/step
#> Epoch 56/100
#> 100/100 - 0s - loss: 0.0046 - 233ms/epoch - 2ms/step
#> Epoch 57/100
#> 100/100 - 0s - loss: 0.0049 - 233ms/epoch - 2ms/step
#> Epoch 58/100
#> 100/100 - 0s - loss: 0.0043 - 230ms/epoch - 2ms/step
#> Epoch 59/100
#> 100/100 - 0s - loss: 0.0055 - 233ms/epoch - 2ms/step
#> Epoch 60/100
#> 100/100 - 0s - loss: 0.0057 - 245ms/epoch - 2ms/step
#> Epoch 61/100
#> 100/100 - 0s - loss: 0.0049 - 232ms/epoch - 2ms/step
#> Epoch 62/100
#> 100/100 - 0s - loss: 0.0048 - 234ms/epoch - 2ms/step
#> Epoch 63/100
#> 100/100 - 0s - loss: 0.0058 - 229ms/epoch - 2ms/step
#> Epoch 64/100
#> 100/100 - 0s - loss: 0.0044 - 218ms/epoch - 2ms/step
#> Epoch 65/100
#> 100/100 - 0s - loss: 0.0044 - 253ms/epoch - 3ms/step
#> Epoch 66/100
#> 100/100 - 0s - loss: 0.0040 - 216ms/epoch - 2ms/step
#> Epoch 67/100
#> 100/100 - 0s - loss: 0.0054 - 226ms/epoch - 2ms/step
#> Epoch 68/100
#> 100/100 - 0s - loss: 0.0043 - 219ms/epoch - 2ms/step
#> Epoch 69/100
#> 100/100 - 0s - loss: 0.0051 - 233ms/epoch - 2ms/step
#> Epoch 70/100
#> 100/100 - 0s - loss: 0.0053 - 216ms/epoch - 2ms/step
#> Epoch 71/100
#> 100/100 - 0s - loss: 0.0041 - 217ms/epoch - 2ms/step
#> Epoch 72/100
#> 100/100 - 0s - loss: 0.0054 - 216ms/epoch - 2ms/step
#> Epoch 73/100
#> 100/100 - 0s - loss: 0.0046 - 234ms/epoch - 2ms/step
#> Epoch 74/100
#> 100/100 - 0s - loss: 0.0044 - 235ms/epoch - 2ms/step
#> Epoch 75/100
#> 100/100 - 0s - loss: 0.0035 - 217ms/epoch - 2ms/step
#> Epoch 76/100
#> 100/100 - 0s - loss: 0.0042 - 216ms/epoch - 2ms/step
#> Epoch 77/100
#> 100/100 - 0s - loss: 0.0049 - 218ms/epoch - 2ms/step
#> Epoch 78/100
#> 100/100 - 0s - loss: 0.0046 - 262ms/epoch - 3ms/step
#> Epoch 79/100
#> 100/100 - 0s - loss: 0.0042 - 241ms/epoch - 2ms/step
#> Epoch 80/100
#> 100/100 - 0s - loss: 0.0038 - 233ms/epoch - 2ms/step
#> Epoch 81/100
#> 100/100 - 0s - loss: 0.0034 - 220ms/epoch - 2ms/step
#> Epoch 82/100
#> 100/100 - 0s - loss: 0.0042 - 234ms/epoch - 2ms/step
#> Epoch 83/100
#> 100/100 - 0s - loss: 0.0046 - 226ms/epoch - 2ms/step
#> Epoch 84/100
#> 100/100 - 0s - loss: 0.0038 - 232ms/epoch - 2ms/step
#> Epoch 85/100
#> 100/100 - 0s - loss: 0.0037 - 230ms/epoch - 2ms/step
#> Epoch 86/100
#> 100/100 - 0s - loss: 0.0042 - 218ms/epoch - 2ms/step
#> Epoch 87/100
#> 100/100 - 0s - loss: 0.0042 - 249ms/epoch - 2ms/step
#> Epoch 88/100
#> 100/100 - 0s - loss: 0.0042 - 232ms/epoch - 2ms/step
#> Epoch 89/100
#> 100/100 - 0s - loss: 0.0030 - 233ms/epoch - 2ms/step
#> Epoch 90/100
#> 100/100 - 0s - loss: 0.0038 - 218ms/epoch - 2ms/step
#> Epoch 91/100
#> 100/100 - 0s - loss: 0.0044 - 234ms/epoch - 2ms/step
#> Epoch 92/100
#> 100/100 - 0s - loss: 0.0030 - 217ms/epoch - 2ms/step
#> Epoch 93/100
#> 100/100 - 0s - loss: 0.0048 - 217ms/epoch - 2ms/step
#> Epoch 94/100
#> 100/100 - 0s - loss: 0.0031 - 219ms/epoch - 2ms/step
#> Epoch 95/100
#> 100/100 - 0s - loss: 0.0042 - 233ms/epoch - 2ms/step
#> Epoch 96/100
#> 100/100 - 0s - loss: 0.0026 - 233ms/epoch - 2ms/step
#> Epoch 97/100
#> 100/100 - 0s - loss: 0.0035 - 218ms/epoch - 2ms/step
#> Epoch 98/100
#> 100/100 - 0s - loss: 0.0040 - 234ms/epoch - 2ms/step
#> Epoch 99/100
#> 100/100 - 0s - loss: 0.0030 - 217ms/epoch - 2ms/step
#> Epoch 100/100
#> 100/100 - 0s - loss: 0.0041 - 217ms/epoch - 2ms/step

Unsurprisingly, Adam does well optimizing a simple neural network such as this.

y_model <- model %>% predict(xdat)
#> 
1/4 [======>.......................] - ETA: 0s
4/4 [==============================] - 0s 1ms/step
#> 
4/4 [==============================] - 0s 1ms/step

plot(xdat,ydat)
lines(xdat,y_model)

SSSoR Attempt

Let us now use SSSoR to optimize the same neural network. First let’s get the layer dimensions of the keras model.

# creation function to assign to weights to keras
weight_dims <- lapply(get_weights(model),dim)
weight_num <- sum(unlist(lapply(weight_dims,prod)))

lowlim=rep(-1,weight_num)
highlim=rep(1,weight_num)

Now we will make a simple function to transform a 1D array into the required Keras weights.

assign_keras_weights <- function(new_array,weight_dims){
  #initialize weight matrix
  weight_mat_new <- list()
  arr_start <- 0
  #Loop through weight matrices
  for (iw in 1:length(weight_dims)){
    matdim = weight_dims[[iw]]
    arr_end = arr_start+prod(matdim)
    weight_mat_new[[iw]] = array(new_array[(arr_start+1):arr_end],dim=matdim)
    arr_start = arr_end
  }
  
  return(weight_mat_new)
}

Now we can establish our model functions and constraints as before.

model_function <-function(params,ls){
  new_weights <- assign_keras_weights(params,weight_dims)
  newmodel <- model
  
  set_weights(newmodel,new_weights)
  
  modval <- newmodel %>% predict(xdat,verbose = F)
  
  return(modval)
}

#Equality Constraints
eq <- function(params){
  constraints = c(
  )
  return(constraints)
}

#Inequality Constraints
# of form <= 0
ineq <- function(params){
  constraints = c(
  )
  return(constraints)
}

Swarm Optimization

Using a Scedastic Surrogate Swarm optimization:

swarm_state <- initialize_swarm(desired_values = ydat,
                                param_len = weight_num,
                                lowlim=lowlim,
                                highlim=highlim,
                                ineq_w = c(),
                                eq_w = c(),
                                config = swarm.config(swarm_size = 100))

for (itt in 1:50){
  tempcontrol = swarm.control(poly_w=0.2,stoch_w = 0.01)
  swarm_state <- step_swarm(desired_values = ydat,
                            swarm_state,
                            ineq_w = c(),
                            eq_w=c(),
                            config = swarm.config(swarm_size=dim(swarm_state$x.p)[1],
                                                  num_cluster=1,
                                                  swarm.control = tempcontrol))
  par(oma=c(0, 0, 0, 5))
  matplot(xdat,swarm_state$pout,type="l",col="grey",xlab = NULL,ylab="value",ylim = c(0,1.2))
  points(xdat,ydat)
  surrogate_out <- surrogate_model(as.vector(swarm_state$poly_rec),swarm_state$polyouts,swarm_state$centersaves,deg=1)
  lines(xdat,surrogate_out,col="blue")
  lines(xdat,model_function(as.vector(swarm_state$poly_rec)),col="red")
  title(paste0("Scedastic Surrogate Swarm Optimization, step = ",itt))
  legend(par('usr')[2], par('usr')[4], bty='n',
         xpd=NA,c("Observed","Modeled","Predicted","Actual"),fill=c("black","grey","blue","red"))
}


print(paste("Number of function evals: ",dim(swarm_state$allxp)[1]))
#> [1] "Number of function evals:  5000"

Search Optimization

We can also use the Scedastic Surrogate Swarm package to do non-swarm searching using surrogates:

search_state <- initialize_search(param_len = weight_num,
                                  lowlim=lowlim,
                                  highlim=highlim,
                                  config = search.config(deg=1))

genstate = F
for (itt in 1:7){
  tempconfig = search.config(gen = genstate,deg=1,
                             search_samples = weight_num*2,
                             search_mag=0.1,
                             revert_best = T)
  search_state <- step_search(desired_values = ydat,search_state,ineq_w = c(),eq_w=c(),config = tempconfig)

  modout = model_function(as.vector(search_state$poly_recs[1,]))
  surrogate_out <- surrogate_model(as.vector(search_state$poly_recs[1,]),search_state$polyouts,search_state$centersaves,deg=1)
  if (mean(abs((modout-surrogate_out)/modout))>10){
    genstate = T
  }else{
    genstate = F
  }
  par(oma=c(0, 0, 0, 5))
  matplot(xdat,search_state$current_pout,type="l",col="grey",xlab = NULL,ylab="value",ylim = c(0,1.2))
  points(xdat,ydat)
  lines(xdat,surrogate_out,col="blue")
  lines(xdat,model_function(as.vector(search_state$poly_recs[1,])),col="red")
  title(paste0("Surrogate Search Optimization, step = ",itt))
  legend(par('usr')[2], par('usr')[4], bty='n', xpd=NA,c("Observed","Modeled","Predicted","Actual"),fill=c("black","grey","blue","red"))
}

print(paste("Number of function evals: ",dim(search_state$xpins)[1]))
#> [1] "Number of function evals:  847"

Search-Initialized Scedastic Surrogate Swarm Optimization

We can then convert these search states and the learned surrogates into the initialization of a Scedastic Surrogate Swarm

#usevec = c(sample(1:dim(search_state$xpins)[1],weight_num),dim(search_state$xpins)[1])
usevec = (dim(search_state$xpins)[1]-100):dim(search_state$xpins)[1]

#Swap into particle swarm
swarm_state <- convert_search_to_swarm(search_state = search_state,
                                       usevec = usevec,
                                       desired_values = ydat,
                                       param_len = weight_num,
                                       lowlim=lowlim,
                                       highlim=highlim,
                                       ineq_w = c(),
                                       eq_w=c(),
                                       config = swarm.config(
                                         swarm_size=length(usevec),
                                         num_cluster=1,
                                         deg=1)
                                       )

for (itt in 1:20){
  tempcontrol = swarm.control(poly_w=0.2,stoch_w = 0.01)
  swarm_state <- step_swarm(desired_values = ydat,
                            swarm_state,
                            ineq_w = c(),
                            eq_w=c(),config = swarm.config(swarm_size=length(usevec),
                                                           num_cluster=1,
                                                           swarm.control = tempcontrol))
  par(oma=c(0, 0, 0, 5))
  matplot(xdat,swarm_state$pout,type="l",col="grey",xlab = NULL,ylab="value",ylim = c(0,1.2))
  points(xdat,ydat)
  surrogate_out <- surrogate_model(as.vector(swarm_state$poly_rec),swarm_state$polyouts,swarm_state$centersaves,deg=1)
  lines(xdat,surrogate_out,col="blue")
  lines(xdat,model_function(as.vector(swarm_state$poly_rec),ls),col="red")
  title(paste0("Search-Initialized Scedastic Surrogate Swarm Optimization, step = ",itt))
  legend(par('usr')[2], par('usr')[4], bty='n',
         xpd=NA,c("Observed","Modeled","Predicted","Actual"),fill=c("black","grey","blue","red"))
}